{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d4d9fdc4",
   "metadata": {},
   "source": [
    "# TLE Fitting\n",
    "\n",
    "Fit a TLE to a set of satellite states (position and velocity)\n",
    "\n",
    "This uses Levenberg-Marquart non-linear least-squares fitting to tune TLE parameters to minimize the difference between the state positions and the SGP4-computed positions\n",
    "\n",
    "Note that a TLE represents states in the TEME frame.  Inputs are rotated into TEME frame from GCRF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd3977b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import satkit as sk\n",
    "import numpy as np\n",
    "import math as m\n",
    "\n",
    "# Create a high-precision state\n",
    "# Altitude for circular orbit\n",
    "altitude = 450e3\n",
    "\n",
    "# Radius & velocity\n",
    "r0 = altitude + sk.consts.earth_radius\n",
    "v0 = m.sqrt(sk.consts.mu_earth / r0)\n",
    "\n",
    "# Inclination\n",
    "inclination = 15 * m.pi / 180.0\n",
    "\n",
    "# Create the state (3D position in meters, 3D velocity in meters / second)\n",
    "state0 = np.array([r0, 0, 0, 0, v0 * m.cos(inclination), v0 * m.sin(inclination)])\n",
    "# Make up an epoch\n",
    "time0 = sk.time(2024, 3, 15, 13, 0, 0)\n",
    "\n",
    "# Propagate the state forward by a day with high-precision propagator\n",
    "res = sk.propagate(state0, time0, time0 + sk.duration(days=1.0))\n",
    "\n",
    "# Get interpolated states every 10 minutes\n",
    "times = [time0 + sk.duration(minutes=i) for i in range(0, 1440, 10)]\n",
    "states = [res.interp(t) for t in times]\n",
    "\n",
    "# Fit the TLE\n",
    "(tle, fitresults) = sk.TLE.fit_from_states(states, times, time0 + sk.duration(days=0.5)) # type: ignore\n",
    "\n",
    "# Print the result\n",
    "print(tle)\n",
    "print(fitresults['success'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b002aeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute position errors (differences between TLE & state)\n",
    "\n",
    "# Get the positions from sgp4\n",
    "(pteme, vteme) = sk.sgp4(tle, times)\n",
    "# Rotate positions from TEME to GCRF frame\n",
    "pgcrf = [sk.frametransform.qteme2gcrf(t) * p for t, p in zip(times, pteme)]\n",
    "# Take difference between state vector and SGP4 positions, and compute norm\n",
    "pdiff = [p - s[0:3] for p, s in zip(pgcrf, states)]\n",
    "pdiff = np.array([np.linalg.norm(p) for p in pdiff])\n",
    "\n",
    "\n",
    "# Plot position errors\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(x=[t.datetime() for t in times], y=pdiff, mode='lines', name='Position Error',\n",
    "                         line=dict(color='black', width=2)))\n",
    "fig.update_layout(title='TLE Fitting Position Errors',\n",
    "                  xaxis_title='Time',\n",
    "                  yaxis_title='Position Error (m)')\n",
    "fig.update_xaxes(showline=True, linewidth=2, linecolor=\"black\", mirror=True)\n",
    "fig.update_yaxes(showline=True, linewidth=2, linecolor=\"black\", mirror=True)\n",
    "fig.update_layout(\n",
    "    xaxis=dict(\n",
    "        gridcolor=\"#dddddd\",\n",
    "        gridwidth=1,\n",
    "    ),\n",
    "    yaxis=dict(\n",
    "        gridcolor=\"#dddddd\",\n",
    "        gridwidth=1,\n",
    "    ),\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}